!>\file cires_ugwpv1_initialize.F90
!!

!===============================
! cu-cires ugwp-scheme
!  initialization of selected 
!  init gw-solvers (1,2,3,4)
!  init gw-source specifications
!  init gw-background dissipation
!==============================
!
! Part-0   specifications of common constants, limiters and "criiical" values
! 
!
    
    module ugwp_common
!
     use machine,  only : kind_phys
			 
     implicit none
     
      real(kind=kind_phys)   ::  pi, pi2, pih, rad_to_deg, deg_to_rad
      real(kind=kind_phys)   ::  arad, p0s     
      real(kind=kind_phys)   ::  grav, grav2, rgrav, rgrav2
      real(kind=kind_phys)   ::  cpd,  rd, rv, fv            
      real(kind=kind_phys)   ::  rdi,  rcpd, rcpd2  
      
      real(kind=kind_phys)   ::  gor,  gr2,  grcp, gocp, rcpdl, grav2cpd
      real(kind=kind_phys)   ::  bnv2min, bnv2max
      real(kind=kind_phys)   ::  dw2min,  velmin, minvel        
      real(kind=kind_phys)   ::  omega1, omega2,   omega3   
      real(kind=kind_phys)   ::  hpscale, rhp, rhp2, rh4, rhp4, khp, hpskm
      real(kind=kind_phys)   ::  mkzmin, mkz2min,  mkzmax, mkz2max, cdmin    
      real(kind=kind_phys)   ::  rcpdt      
           
!      real(kind=kind_phys), parameter ::  grav2 = grav + grav
!      real(kind=kind_phys), parameter ::  rgrav = 1.0/grav, rgrav2= rgrav*rgrav           
!      real(kind=kind_phys), parameter ::  rdi  = 1.0 / rd, rcpd = 1./cpd, rcpd2 = 0.5/cpd  
!      real(kind=kind_phys), parameter ::  gor  = grav/rd, rcpdt = 1./(cp*dtp)
      
!      real(kind=kind_phys), parameter ::  gr2  = grav*gor
!      real(kind=kind_phys), parameter ::  grcp = grav*rcpd, gocp = grcp
!      real(kind=kind_phys), parameter ::  rcpdl    = cpd*rgrav                         ! 1/[g/cp]  == cp/g
!      real(kind=kind_phys), parameter ::  grav2cpd = grav*grcp                         !  g*(g/cp)= g^2/cp
!      real(kind=kind_phys), parameter ::  pi2 = 2.*pi, pih = .5*pi  
!      real(kind=kind_phys), parameter ::  rad_to_deg=180.0/pi, deg_to_rad=pi/180.0
!               
!      real(kind=kind_phys), parameter ::  bnv2min = (pi2/1800.)*(pi2/1800.)
!      real(kind=kind_phys), parameter ::  bnv2max = (pi2/30.)*(pi2/30.)      
!      real(kind=kind_phys), parameter :: dw2min=1.0,  velmin=sqrt(dw2min), minvel = 0.5     
!      real(kind=kind_phys), parameter :: omega1 = pi2/86400., omega2 = 2.*omega1, omega3 = 3.*omega1     
!   
!      real(kind=kind_phys), parameter :: hpscale= 7000., rhp=1./hpscale, rhp2=.5*rhp, rh4 = 0.25*rhp 
!      real(kind=kind_phys), parameter :: mkzmin = pi2/80.0e3, mkz2min = mkzmin*mkzmin
!      real(kind=kind_phys), parameter :: mkzmax = pi2/500.,   mkz2max = mkzmax*mkzmax 
!      real(kind=kind_phys), parameter :: cdmin  = 2.e-2/mkzmax 
!      real(kind=kind_phys), parameter ::  pi = 4.*atan(1.0),
!      real(kind=kind_phys), parameter ::  grav =9.81, cpd = 1004.
!      real(kind=kind_phys), parameter ::  rd = 287.0 , rv =461.5  
!      real(kind=kind_phys), parameter ::  fv   = rv/rd - 1.0    
!      real(kind=kind_phys), parameter ::  arad = 6370.e3
        
     end module ugwp_common
     
      subroutine init_nazdir(naz,  xaz,  yaz)
      
      use machine,     only : kind_phys
      use ugwp_common, only :  pi2
      
      implicit none
      
      integer :: naz
      real(kind=kind_phys), dimension(naz) :: xaz,  yaz
      integer :: idir
      real(kind=kind_phys)    :: phic, drad
      
      drad  = pi2/float(naz)
      if (naz.ne.4) then     
        do idir =1, naz
         Phic = drad*(float(idir)-1.0)
         xaz(idir) = cos(Phic)
         yaz(idir) = sin(Phic)
        enddo
      else 
!                                  if (naz.eq.4) then
          xaz(1) = 1.0     !E
          yaz(1) = 0.0
          xaz(2) = 0.0     
          yaz(2) = 1.0     !N
          xaz(3) =-1.0     !W
          yaz(3) = 0.0
          xaz(4) = 0.0
          yaz(4) =-1.0     !S
      endif      
      end  subroutine init_nazdir     
!
!
!===================================================
!
!Part-1 init =>   wave dissipation + RFriction
!
!===================================================
     subroutine init_global_gwdis(levs, zkm, pmb, kvg, ktg, krad, kion, me, master)
!				      
     use machine ,      only : kind_phys 
     use ugwp_common,   only : pih, pi  
     
     implicit none
     integer , intent(in)                 :: me, master
     integer , intent(in)                 :: levs
     real(kind=kind_phys), intent(in)                     :: zkm(levs), pmb(levs)    ! in km-Pa
     real(kind=kind_phys), intent(out), dimension(levs+1) :: kvg, ktg, krad, kion
!
!locals + data
!
     integer :: k
     real(kind=kind_phys), parameter :: vusurf  = 2.e-5
     real(kind=kind_phys), parameter :: musurf  = vusurf/1.95
     real(kind=kind_phys), parameter :: hpmol   = 7.0
!
     real(kind=kind_phys), parameter :: kzmin   =  0.1
     real(kind=kind_phys), parameter :: kturbo  = 100.
     real(kind=kind_phys), parameter :: zturbo  = 130.
     real(kind=kind_phys), parameter :: zturw   = 30.
     real(kind=kind_phys), parameter :: inv_pra = 3.  !kt/kv =inv_pr
!
     real(kind=kind_phys), parameter :: alpha  = 1./86400./15.   ! height variable see Zhu-1993 from 60-days => 6 days
     real(kind=kind_phys)  ::  pa_alp  = 750.                    ! super-RF parameters from FV3-dycore GFSv17/16 sett
     real(kind=kind_phys)  :: tau_alp  = 10.                     ! days   (750 Pa /10days)
!     
     real(kind=kind_phys), parameter :: kdrag  = 1./86400./30.   !parametrization for WAM ion drag as e-density function
     real(kind=kind_phys), parameter :: zdrag  = 100.
     real(kind=kind_phys), parameter :: zgrow  = 50.
!
     real(kind=kind_phys) ::   vumol, mumol, keddy, ion_drag
     real(kind=kind_phys) ::   rf_fv3, rtau_fv3, ptop, pih_dlog       
!
     real(kind=kind_phys) ::  ae1 ,ae2
!     
      
     ptop = pmb(levs)  
     rtau_fv3 = 1./86400./tau_alp
     pih_dlog = pih/log(pa_alp/ptop)

     do k=1, levs
      ae1 =    zkm(k)/hpmol
      vumol    = vusurf*exp(ae1)
      mumol    = musurf*exp(ae1)
      ae2 = -((zkm(k)-zturbo) /zturw)**2
      keddy    = kturbo*exp(ae2)

      kvg(k)   = vumol + keddy
      ktg(k)   = mumol + keddy*inv_pra

      krad(k)  = alpha
!
      ion_drag = kdrag
!
      kion(k)  = ion_drag!      
! add  Rayleigh_Super of FV3  for pmb < pa_alp
!
      if (pmb(k) .le. pa_alp) then
       rf_fv3=rtau_fv3*sin(pih_dlog*log(pa_alp/pmb(k)))**2
       krad(k) = krad(k) + rf_fv3 
       kion(k) = kion(k) + rf_fv3   

      endif

!      write(6,132) zkm(k), kvg(k), kvg(k)*(6.28/5000.)**2,  kion(k) 
     enddo

      k= levs+1
      kion(k) = kion(k-1)
      krad(k) = krad(k-1)
      kvg(k)  =  kvg(k-1)
      ktg(k)  =  ktg(k-1)
      
!      if (me == master) then
!	  write(6, * ) '  zkm(k), kvg(k), kvg(k)*(6.28/5000.)**2,  kion(k) ... init_global_gwdis'     
!        do k=1, levs, 1
!	  write(6,132) zkm(k), kvg(k), kvg(k)*(6.28/5000.)**2,  kion(k), pmb(k) 
!	enddo      
!      endif
!
! 132  format( 2x, F8.3,' dis-scales:', 4(2x, E10.3))    
                                                          
     end subroutine init_global_gwdis
!
! ========================================================================
! Part 2 - sources
!      wave  sources
! ========================================================================
!
!    ugwp_oro_init
!
!=========================================================================
     module ugwp_oro_init
     use machine ,    only : kind_phys
     use ugwp_common, only : bnv2min, grav, grcp, fv, grav, cpd, grcp, pi
     use ugwp_common, only : mkzmin, mkz2min
     implicit none
!  
! constants and "crirtical" values to run oro-mtb_gw physics
!
!
      real(kind=kind_phys),      parameter :: hncrit=9000.   ! max value in meters for elvmax
      real(kind=kind_phys),      parameter :: hminmt=50.     ! min mtn height (*j*)
      real(kind=kind_phys),      parameter :: sigfac=4.0     ! mb3a expt test for elvmax factor
      real(kind=kind_phys),      parameter :: hpmax=2500.0
      real(kind=kind_phys),      parameter :: hpmin=25.0      
!
!
      real(kind=kind_phys),      parameter :: minwnd=1.0     ! min wind component (*j*)
      real(kind=kind_phys),      parameter :: dpmin=5000.0   ! minimum thickness of the reference layer in pa

      
      character(len=8)  ::  strver  = 'gfs_2018'
      character(len=8)  ::  strbase = 'gfs_2018'
      
      real(kind=kind_phys), parameter   ::  rimin=-10., ric=0.25
      
      real(kind=kind_phys), parameter   ::  frmax=10., frc =1.0, frmin =0.01
      real(kind=kind_phys), parameter   ::  ce=0.8,   ceofrc=ce/frc, cg=0.5
      real(kind=kind_phys), parameter   ::  gmax=1.0, veleps=1.0, factop=0.5!
      real(kind=kind_phys), parameter   ::  efmin=0.5,    efmax=10.0
      
      real(kind=kind_phys), parameter   ::  rlolev=50000.0
      integer,              parameter   ::  mdir = 8
      real(kind=kind_phys), parameter   ::  fdir=mdir/(8.*atan(1.0))
      real(kind=kind_phys), parameter   ::  zpgeo=2.*atan(1.0)
      
      integer nwdir(mdir)
      data nwdir/6,7,5,8,2,3,1,4/
      save nwdir
      
      real(kind=kind_phys), parameter   :: odmin = 0.1, odmax = 10.0
      real(kind=kind_phys), parameter   :: fcrit_sm  = 0.7, fcrit_sm2  = fcrit_sm * fcrit_sm
      real(kind=kind_phys), parameter   :: fcrit_gfs = 0.7, fcrit_v1 = 0.7
      real(kind=kind_phys), parameter   :: fcrit_mtb = 0.7

      real(kind=kind_phys),  parameter  :: zbr_pi  = zpgeo
      real(kind=kind_phys),  parameter  :: zbr_ifs = zpgeo
 
!  
 
      real(kind=kind_phys),   parameter ::   kxoro=6.28e-3/200.    !
      real(kind=kind_phys),   parameter ::   coro = 0.0
      integer,parameter ::    nridge=2
      real(kind=kind_phys),   parameter ::   sigma_std=1./100., gamm_std=1.0       
 
      real(kind=kind_phys)    ::  cdmb                      ! scale factors for mtb
      real(kind=kind_phys)    ::  cleff                     ! scale factors for orogw
      
      integer ::  nworo                     ! number of waves
      integer ::  nazoro                    ! number of azimuths
      integer ::  nstoro                    ! flag for stochastic launch above SG-peak

      
!------------------------------------------------------------------------------
!    small-scale orography parameters  for TOFD of Beljaars et al., 2004, QJRMS
!    SA-option can be controlled by  Integral limits of fluxes
!    in B2004: klow = 0.003 1/m ~ 2km and kinf ~ 6.28/10/(Z1)~< 1 km => meters
!    these limits can change strength of TOFD... choice of k0tr ~1/10 km (10km ~dx of C768)
!    kmax = kdis_pbl
!------------------------------------------------------------------------------
     real(kind=kind_phys), parameter     :: kmax = 6.28/(10.*25.)           ! max k-tofd
     real(kind=kind_phys), parameter     :: k1tr = 6.28/(2100)              ! max k-transition from -1.9/slope to -2.8/slope
     real(kind=kind_phys), parameter     :: kflt = 6.28/(18.e3)             !
     real(kind=kind_phys), parameter     :: k0tr = 6.28/(10.e3)             ! min k-tofd     
     real(kind=kind_phys), parameter     :: nk1tr = 2.8
     real(kind=kind_phys), parameter     :: nk0tr = 1.9
     real(kind=kind_phys), parameter     :: a1_tofd =  kflt ** nk1tr *1.e3
     real(kind=kind_phys), parameter     :: a2_tofd =  k1tr ** (nk0tr-nk1tr)    
     real(kind=kind_phys), parameter     :: fix_tofd = 2.* 0.005 * 12 *0.6     !value= 0.072
!     
! B2004 scheme is based on the empirical vertical profile of the tofd divergence: 
!       Ax_tofd(Z)=exp(-[Z/ze_tofd]^3/2) / Z^1.2.....
! TOFD-flux/TMS-flux must dissipate due to PBL-diffusion with spectral damping
! Here we can enhance TOFD-impact by selecting k0tr and kmax limits
!      as functions of resolution and PBL-dissipation
!                            
     integer, parameter  :: n_tofd = 2                                      ! depth of SSO for TOFD compared with Zpbl
     real(kind=kind_phys), parameter     :: const_tofd = 0.0759             ! alpha*beta*Cmd*Ccorr*2.109 = 12.*1.*0.005*0.6*2.109 = 0.0759
     real(kind=kind_phys), parameter     :: ze_tofd    = 1500.0             ! BJ's z-decay in meters, 1.5 km
     real(kind=kind_phys), parameter     :: a12_tofd   = 0.0002662*0.005363 ! BJ's k-spect const for sigf2 * a1*a2*exp(-[z/zdec]**1.5]
     real(kind=kind_phys), parameter     :: ztop_tofd  = 3.*ze_tofd         ! no TOFD > this height 4.5 km
!------------------------------------------------------------------------------
!

      contains
!
      subroutine init_oro_gws(nwaves, nazdir, nstoch, effac, &
                              lonr, kxw )
!
!
      integer :: nwaves, nazdir, nstoch
      integer :: lonr			       
      real(kind=kind_phys)    :: cdmbX
      real(kind=kind_phys)    :: kxw
      real(kind=kind_phys)    :: effac        ! it is analog of cdmbgwd(2) for GWs, off for now
!-----------------------------! GFS-setup for cdmb & cleff
! cdmb =  4.0 * (192.0/IMX)  
! cleff = 0.5E-5 / SQRT(IMX/192.0)  = 0.5E-5*SQRT(192./IMX)
!
      real(kind=kind_phys), parameter :: lonr_refmb =  4.0 * 192.0
      real(kind=kind_phys), parameter :: lonr_refgw =  192.0
      real(kind=kind_phys), parameter :: cleff_ref  = 0.5e-5           !       1256 km = 10 * 125 km ???
      
! copy  to "ugwp_oro_init"  =>  nwaves, nazdir, nstoch
 
      nworo  =  nwaves
      nazoro =  nazdir
      nstoro =  nstoch

      cdmbX = lonr_refmb/float(lonr)
      
      cdmb  = cdmbX
      cleff = cleff_ref * sqrt(lonr_refgw/float(lonr))   !* effac
!
      end subroutine init_oro_gws
!

    end module ugwp_oro_init
! =========================================================================
!
!    ugwp_conv_init
!
!=========================================================================
    module ugwp_conv_init
    
     use machine ,              only : kind_phys

     
     implicit none
      real(kind=kind_phys)    ::  eff_con                   ! scale factors for conv GWs
      integer ::  nwcon                     ! number of waves
      integer ::  nazcon                    ! number of azimuths
      integer ::  nstcon                    ! flag for stochastic choice of launch level above Conv-cloud
      real(kind=kind_phys)    ::  con_dlength
      real(kind=kind_phys)    ::  con_cldf

      real(kind=kind_phys), parameter    :: cmin  =  5  !2.5
      real(kind=kind_phys), parameter    :: cmax  = 95. !82.5
      real(kind=kind_phys), parameter    :: cmid  = 22.5
      real(kind=kind_phys), parameter    :: cwid  = cmid
      real(kind=kind_phys), parameter    :: bns   = 2.e-2, bns2 = bns*bns, bns4=bns2*bns2
      real(kind=kind_phys), parameter    :: mstar = 6.28e-3/2.  !  2km
      real(kind=kind_phys)               :: dc

      real(kind=kind_phys), allocatable  :: ch_conv(:),  spf_conv(:)
      real(kind=kind_phys), allocatable  :: xaz_conv(:), yaz_conv(:)
     contains
!
     subroutine init_conv_gws(nwaves, nazdir, nstoch, effac, lonr, kxw)
!
     use ugwp_common,  only : pi2, arad

     implicit none
     
 
      integer :: nwaves, nazdir, nstoch
      integer :: lonr
!      
! ccpp
!      
      
      real(kind=kind_phys)    :: kxw,  effac
      real(kind=kind_phys)    :: work1 = 0.5
      real(kind=kind_phys)    :: chk, tn4, snorm
      integer :: k

      nwcon    = nwaves
      nazcon   = nazdir
      nstcon   = nstoch
      eff_con  = effac

      con_dlength = pi2*arad/float(lonr)      
!
! allocate & define spectra in "selected direction": "dc" "ch(nwaves)"
!
       if (.not. allocated(ch_conv))  allocate (ch_conv(nwaves))
       if (.not. allocated(spf_conv)) allocate (spf_conv(nwaves))
       if (.not. allocated(xaz_conv)) allocate (xaz_conv(nazdir))
       if (.not. allocated(yaz_conv)) allocate (yaz_conv(nazdir))

      dc = (cmax-cmin)/float(nwaves-1)
!
! we may use different spectral "shapes"
! for example FVS-93 "Desabeius"
!  E(s=1, t=3,m, w, k) ~ m^s/(m*^4 + m^4) ~ m^-3 saturated tail
!
      do k = 1,nwaves
         chk         = cmin + (k-1)*dc
         tn4         = (mstar*chk)**4
         ch_conv(k)  =  chk
         spf_conv(k) =  bns4*chk/(bns4+tn4)
      enddo

      snorm = sum(spf_conv)
      spf_conv = spf_conv/snorm*1.5
 
      call init_nazdir(nazdir,  xaz_conv,  yaz_conv)
     end subroutine init_conv_gws


    end module ugwp_conv_init
!=========================================================================
!
!    ugwp_fjet_init
!
!=========================================================================

   module ugwp_fjet_init
   use machine ,              only : kind_phys   

   
   
      implicit none
      real(kind=kind_phys)    ::  eff_fj                     ! scale factors for conv GWs
      integer ::  nwfj                       ! number of waves
      integer ::  nazfj                      ! number of azimuths
      integer ::  nstfj                      ! flag for stochastic choice of launch level above Conv-cloud
!
      real(kind=kind_phys), parameter    ::  fjet_trig=0.    ! if ( abs(frgf) > fjet_trig ) launch GW-packet

      
      real(kind=kind_phys), parameter    :: cmin =  2.5
      real(kind=kind_phys), parameter    :: cmax = 67.5
      real(kind=kind_phys)               :: dc
      real(kind=kind_phys), allocatable  :: ch_fjet(:) , spf_fjet(:)
      real(kind=kind_phys), allocatable  :: xaz_fjet(:), yaz_fjet(:)
     contains
     
     subroutine init_fjet_gws(nwaves, nazdir, nstoch, effac,lonr, kxw) 

     use ugwp_common,  only : pi2, arad			      
			      
     implicit none

      integer :: nwaves, nazdir, nstoch
      integer :: lonr
      real(kind=kind_phys)    :: kxw,  effac , chk
 
      integer :: k

      nwfj   =  nwaves
      nazfj  =  nazdir
      nstfj  =  nstoch
      eff_fj =  effac

       if (.not. allocated(ch_fjet))  allocate (ch_fjet(nwaves))
       if (.not. allocated(spf_fjet)) allocate (spf_fjet(nwaves))
       if (.not. allocated(xaz_fjet)) allocate (xaz_fjet(nazdir))
       if (.not. allocated(yaz_fjet)) allocate (yaz_fjet(nazdir))
 
      dc = (cmax-cmin)/float(nwaves-1)
      do k = 1,nwaves
         chk         = cmin + (k-1)*dc
         ch_fjet(k)  =  chk
         spf_fjet(k) =  1.0
      enddo
      call init_nazdir(nazdir,  xaz_fjet,  yaz_fjet)

     end subroutine init_fjet_gws

    end module ugwp_fjet_init
!
!=========================================================================
!
!
     module ugwp_okw_init
!=========================================================================
       use machine ,              only : kind_phys      
     
       implicit none

      real(kind=kind_phys)    ::  eff_okw                     ! scale factors for conv GWs
      integer ::  nwokw                       ! number of waves
      integer ::  nazokw                      ! number of azimuths
      integer ::  nstokw                      ! flag for stochastic choice of launch level above Conv-cloud
!
      real(kind=kind_phys), parameter    ::  okw_trig=0.      ! if ( abs(okwp) > okw_trig ) launch GW-packet

      real(kind=kind_phys), parameter    :: cmin =  2.5
      real(kind=kind_phys), parameter    :: cmax = 67.5
      real(kind=kind_phys)               :: dc
      real(kind=kind_phys), allocatable  :: ch_okwp(:),   spf_okwp(:)
      real(kind=kind_phys), allocatable  :: xaz_okwp(:),  yaz_okwp(:)

     contains
!

			     
     subroutine init_okw_gws(nwaves, nazdir, nstoch, effac, lonr, kxw)
     use ugwp_common,  only : pi2, arad

     implicit none

      integer :: nwaves, nazdir, nstoch
      integer :: lonr
      real(kind=kind_phys)    :: kxw,  effac , chk
 
      integer :: k

      nwokw   =  nwaves
      nazokw  =  nazdir
      nstokw  =  nstoch
      eff_okw =  effac

       if (.not. allocated(ch_okwp))  allocate (ch_okwp(nwaves))
       if (.not. allocated(spf_okwp)) allocate (spf_okwp(nwaves))
       if (.not. allocated(xaz_okwp)) allocate (xaz_okwp(nazdir))
       if (.not. allocated(yaz_okwp)) allocate (yaz_okwp(nazdir))
      dc = (cmax-cmin)/float(nwaves-1)
      do k = 1,nwaves
         chk =  cmin + (k-1)*dc
         ch_okwp(k) = chk
         spf_okwp(k) = 1.
      enddo
  
      call init_nazdir(nazdir,  xaz_okwp,  yaz_okwp)
!
     end subroutine init_okw_gws
 
     end module ugwp_okw_init

!=============================== end of GW  sources
!
!  init specific  gw-solvers (1,2,3,4)
!
!===============================
!  Part -3  init  wave solvers
!===============================

  module ugwp_lsatdis_init
     use machine ,              only : kind_phys    
     implicit none

      integer  :: nwav, nazd
      integer  :: nst
      real(kind=kind_phys)     :: eff
      integer, parameter  :: incdim = 4, iazdim = 4
!
     contains

     subroutine initsolv_lsatdis(me, master,  nwaves, nazdir, nstoch, effac, do_physb, kxw)

     implicit none
!
     integer  :: me, master
     integer  :: nwaves, nazdir
     integer  :: nstoch
     real(kind=kind_phys)     :: effac
     logical  :: do_physb
     real(kind=kind_phys)     :: kxw
!
!locals: define azimuths and Ch(nwaves) - domain when physics-based soureces
!                                          are not actibve
!
     integer :: inc, jk, jl, iazi, i, j, k

     if( nwaves == 0 .or. nstoch == 1 ) then
!                                redefine from the default
       nwav = incdim
       nazd = iazdim
       nst  = 0
       eff  = 1.0
     else
!                                from input_nml multi-wave spectra
       nwav = nwaves
       nazd = nazdir
       nst  = nstoch
       eff  = effac
     endif
!
     end subroutine initsolv_lsatdis
!
  end module ugwp_lsatdis_init
!
!
  module ugwp_wmsdis_init
  
     use machine ,    only : kind_phys    
     use ugwp_common, only : arad, pi, pi2, hpscale, rhp, rhp2, rh4, omega2 
     use ugwp_common, only : bnv2max,  bnv2min, minvel 
     use ugwp_common, only : mkzmin,  mkz2min, mkzmax, mkz2max, ucrit => cdmin
     
    implicit none

      real(kind=kind_phys),     parameter   :: maxdudt = 250.e-5, maxdtdt=15.e-2
      real(kind=kind_phys),     parameter   :: dked_min =0.01,    dked_max=250.0
   
      real(kind=kind_phys),     parameter   :: gptwo=2.0
 
      real(kind=kind_phys) ,     parameter  :: bnfix  =  6.28/300., bnfix2= bnfix * bnfix
      real(kind=kind_phys) ,     parameter  :: bnfix4 =  bnfix2 * bnfix2  
      real(kind=kind_phys) ,     parameter  :: bnfix3 =  bnfix2 * bnfix        
!
! make parameter list that will be passed to SOLVER
! 
      integer  , parameter  :: iazidim=4       ! number of azimuths
      integer  , parameter  :: incdim=25       ! number of discrete cx - spectral elements in launch spectrum
 
      real(kind=kind_phys) ,     parameter  :: zcimin = 2.5
      real(kind=kind_phys) ,     parameter  :: zcimax = 125.0
      real(kind=kind_phys) ,     parameter  :: zgam   = 0.25
!
! Verical spectra
!
      real(kind=kind_phys) ,     parameter  :: pind_wd = 5./3.
      real(kind=kind_phys) ,     parameter  :: sind_kz = 1.
      real(kind=kind_phys) ,     parameter  :: tind_kz = 3.    
      real(kind=kind_phys) ,     parameter  :: stind_kz = sind_kz + tind_kz  
!
! copies from kmob_ugwp namelist
!      
      real(kind=kind_phys)                  :: nslope              ! the GW sprctral slope at small-m             
      real(kind=kind_phys)                  :: lzstar
      real(kind=kind_phys)                  :: lzmin
      real(kind=kind_phys)                  :: lzmax      
      real(kind=kind_phys)                  :: lhmet
      real(kind=kind_phys)                  :: tamp_mpa            !amplitude for GEOS-5/MERRA-2
      real(kind=kind_phys)                  :: tau_min             ! min of GW MF 0.25 mPa                        
      integer               :: ilaunch
      real(kind=kind_phys)                  :: gw_eff

      real(kind=kind_phys)                  :: v_kxw, rv_kxw, v_kxw2 
     
     
 
!===========================================================================
      integer  :: nwav, nazd, nst
      real(kind=kind_phys)     :: eff
 
      real(kind=kind_phys)                :: zaz_fct, zms
      real(kind=kind_phys), allocatable   :: zci(:), zci4(:), zci3(:),zci2(:), zdci(:)
      real(kind=kind_phys), allocatable   :: zcosang(:), zsinang(:)
      real(kind=kind_phys), allocatable   :: lzmet(:), czmet(:), mkzmet(:), dczmet(:), dmkz(:) 

!
! GW-eddy constants for wave-mode dissipation by background and stability of
!         "final" flow after application of GW-effects
! 
      real(kind=kind_phys), parameter :: iPr_pt = 0.5
      real(kind=kind_phys), parameter :: lturb = 30., sc2  = lturb*lturb     ! stable on 80-km TL lmix ~ 500 met.
      real(kind=kind_phys), parameter :: ulturb=150., sc2u = ulturb* ulturb  ! unstable
      real(kind=kind_phys), parameter :: ric =0.25
      real(kind=kind_phys), parameter :: rimin = -10., prmin = 0.25
      real(kind=kind_phys), parameter :: prmax = 4.0
!
      contains
!============================================================================
     subroutine initsolv_wmsdis(me, master,  nwaves, nazdir, nstoch, effac, do_physb, kxw, version)
 
!        call initsolv_wmsdis(me, master, knob_ugwp_wvspec(2), knob_ugwp_azdir(2), &
!         knob_ugwp_stoch(2), knob_ugwp_effac(2), do_physb_gwsrcs, kxw,version)
!
     implicit none
!
!input-control for solvers:
!      nwaves, nazdir, nstoch, effac, do_physb, kxw
!
!
     integer, intent(in)  :: me, master, nwaves, nazdir, nstoch
     integer, intent(in)  :: version 
        
     real(kind=kind_phys), intent(in)     :: effac, kxw
     logical,  intent(in) :: do_physb
      
!
!locals
!
      real(kind=kind_phys)     :: dlzmet   
      real(kind=kind_phys)     :: cstar,rcstar, nslope3, fnorm, zcin     
            
      integer  :: inc, jk, jl, iazi
!
      real(kind=kind_phys) :: zang, zang1, znorm
      real(kind=kind_phys) :: zx1, zx2, ztx, zdx, zxran, zxmin, zxmax, zx, zpexp
      real(kind=kind_phys) :: fpc, fpc_dc
      real(kind=kind_phys) :: ae1,ae2
     if( nwaves == 0) then
!
!     redefine from the deafault
!
       nwav   = incdim
       nazd   = iazidim
       nst    = 0
       eff    = 1.0      
       gw_eff = eff
     else
!
! from input.nml
!
       nwav   = nwaves
       nazd   = nazdir
       nst    = nstoch
       gw_eff = effac
     endif  

      
      v_kxw  = kxw ; v_kxw2 = v_kxw*v_kxw 
      rv_kxw = 1./v_kxw

      allocate ( zci(nwav),  zci4(nwav), zci3(nwav),zci2(nwav), zdci(nwav)  )
      allocate ( zcosang(nazd), zsinang(nazd) )
      allocate (lzmet(nwav), czmet(nwav), mkzmet(nwav), dczmet(nwav), dmkz(nwav) )   

!      if (me == master) then
!         print *, 'ugwp_v1/v0: init_gw_wmsdis_control '
!  
!         print *, 'ugwp_v1/v0: WMS_DIS launch layer ',    ilaunch
!         print *, 'ugwp_v1/v0: WMS_DIS tot_mflux in mpa', tamp_mpa*1000.
!         print *, 'ugwp_v1/v0: WMS_DIS lhmet in km  ' ,   lhmet*1.e-3      	 
!       endif

       zpexp = gptwo * 0.5                    ! gptwo=2 , zpexp = 1.

!
!      set up azimuth directions and some trig factors
!
!
       zang = pi2 / float(nazd)

! get normalization factor to ensure that the same amount of momentum
! flux is directed (n,s,e,w) no mater how many azimuths are selected.
!
       znorm = 0.0
       do iazi=1, nazd
         zang1         = (iazi-1)*zang
         zcosang(iazi) = cos(zang1)
         zsinang(iazi) = sin(zang1)
         znorm         = znorm + abs(zcosang(iazi))
       enddo
!      zaz_fct = 1.0
       zaz_fct = 2.0 / znorm            ! correction factor for azimuthal sums

!       define coordinate transform for "Ch"   ....x = 1/c stretching transform
!       ----------------------------------------------- 
! 
! x=1/Cphase transform
! Scinocca 2003.   x = 1/c stretching transform
!
          zxmax = 1.0 / zcimin
          zxmin = 1.0 / zcimax
          zxran = zxmax - zxmin
          zdx   = zxran / real(nwav-1)                            ! dkz
!
          ae1=zxran/zgam
          zx1   = zxran/(exp(ae1)-1.0 )                    ! zgam =1./4.
          zx2   = zxmin - zx1

!
! computations for zci =1/zx, stretching "accuracy" is not "accurate" spectra transform
!  it represents additional "empirical" redistribution of "spectral" mode in C-space
!                               
         zms = pi2 / lzstar

        do inc=1, nwav
          ztx = real(inc-1)*zdx+zxmin
	  ae1 = (ztx-zxmin)/zgam
          zx  = zx1*exp(ae1)+zx2                            !eq.(29-30),Scinocca-2003
          zci(inc)  = 1.0 /zx                                            !
          zdci(inc) = zci(inc)**2*(zx1/zgam)*exp(ae1)*zdx   !
          zci4(inc) = (zms*zci(inc))**4
          zci2(inc) = (zms*zci(inc))**2
          zci3(inc) = (zms*zci(inc))**3
        enddo
!
!
!  alternatuve lzmax-lzmin without  x=1/c transform
!
!
        if (version == 1) then 
	
	dlzmet = (lzmax-lzmin)/ real(nwav-1)    
         do inc=1, nwav
          lzmet(inc) = lzmin + (inc-1)*dlzmet
	  mkzmet(inc) = pi2/lzmet(inc)
          zci(inc) =lzmet(inc)/(pi2/bnfix)
          zci4(inc) = (zms*zci(inc))**4
          zci2(inc) = (zms*zci(inc))**2
          zci3(inc) = (zms*zci(inc))**3	
	  
	 enddo

	  zdx = (zci(nwav)-zci(1))/ real(nwav-1)  
          do inc=1, nwav	    	  
          zdci(inc) = zdx     
	  enddo   
	 
	  cstar = bnfix/zms
	  rcstar = 1./cstar	
	ENDIF                               !   if (version == 1) then 	 
	 
	RETURN	  	 
!===================  Diag prints after return  ====================
         if (me == master) then
           print *
           print *,  'ugwp_v0: zcimin=' , zcimin
           print *,  'ugwp_v0: zcimax=' , zcimax
           print *,  'ugwp_v0: zgam=   ', zgam                
           print *
           
!          print *, ' ugwp_v1  nslope=', nslope
           print *
           print *,  'ugwp_v1: zcimin/zci=' , maxval(zci)
           print *,  'ugwp_v1: zcimax/zci=' , minval(zci)
           print *,  'ugwp_v1: cd_crit=', ucrit            
           print *,  'ugwp_v1: launch_level',  ilaunch
           print *, ' ugwp_v1  lzstar=', lzstar
           print *, ' ugwp_v1  nslope=', nslope

           print *
	   nslope3=nslope+3.0
         do inc=1, nwav	    	  
           zcin =zci(inc)*rcstar
	   fpc =  rcstar*(zcin*zcin)/(1.+ zcin**nslope3)	   	    
           fpc_dc = fpc * zdci(inc)
	   write(6,111)  inc, zci(inc), zdci(inc),ucrit, fpc, fpc_dc, 6.28e-3/bnfix*zci(inc)
          enddo
         endif
	 

	 
 111     format( 'wms-zci', i4, 7 (3x, F8.3))

     end subroutine initsolv_wmsdis
!
!
  end module ugwp_wmsdis_init
